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The mechanism of Mott transition in the Hubbard model on the square lattice is studied 
without explicit introduction of magnetic and superconducting correlations, using a variational 
Monte Carlo method. In the trial wave functions, we consider various types of binding factors 
between a doubly-occupied site (doublon, D) and an empty site (holon, H), like a long-range 
type as well as a conventional nearest-neighbor type, and add independent long-range D-D 
(H-H) factors. It is found that a wide choice of D-H binding factor leads to Mott transitions at 
critical values near the band width. We renew the D-H binding picture of Mott transitions by 
introducing two characteristic length scales, the D-H binding length Ion and the minimum D-D 
distance <dd, which we appropriately estimate. A Mott transition takes place at Ion ~ Idd- In 
the metallic regime (^dh > ^dd), the domains of D-H pairs overlap with one another, thereby 
doublons and holons can move independently by exchanging the partners one after another. In 
contrast, the D-D factors give only a minor contribution to the Mott transition. 

KEYWORDS: Mott transition, Hubbard model, variational Monte Carlo, doublon-holon binding, square 
lattice 



1. Introduction 

The Mott transitions 1 - 1 free from the spin degree of 
freedom were recently realized in ultracold bosonic atoms 
on optical lattices. 2 ~ 5 ' In the Bose Hubbard models, 6 ' 
which faithfully describe these systems, 7 ' 8 ' superfluid- 
insulator (Mott) transitions occur when kinetic and in- 
teraction energies are competitive. 9 ~ 16 ' It follows that 
the essence of "Mott physics" can be separated from 
magnetic metal-insulator transitions, 17 ' often arising in 
weakly correlated regimes in half-filled-band electronic 
systems, especially, with good nesting conditions. Aside 
from the spinless bosons, systems like organic super- 
conductors k-(BEDT-TTF) salts 18 ' and the cuprate su- 
perconductors as doped Mott insulators 19-21 ' require a 
deep understanding about the mechanism of nonmag- 
netic Mott transitions. Thus, it is significant to shed 
light on the Mott transition between virtual paramag- 
netic phases in the Hubbard model on the square lattice, 
which, though, actually has an antiferromagnetic (AF) 
long-range order for any finite value of positive U/t (U: 
onsite correlation strength; t: hopping integral). 22 ' 23 ' 

The variation theory 24 ' 25 ' has long been one of the 
main streams to study ground-state properties of the 
Mott transition. In particular, the variational Monte 
Carlo (VMC) approaches 26-28 ' are effective for its re- 
liability in dealing with the local correlation and wide 
applicability. In the previous VMC studies related to 
paramagnetic Mott transitions, the following properties 
have been clarified, (i) The well-known Gutzwiller wave 
function, 24 ' with only onsite correlation, does not un- 
dergo a Mott transition for finite U/t and in finite di- 
mensions. Namely, the Brinkman-Rice transition 25 ' 29 ' is 
unreal. 28,30 ' (ii) Wave functions with short-range inter- 
site attractive correlations between a doubly-occupied 
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site (doublon, D) and an empty site (holon, H), 31 ~ 33 ' 
which are minus and plus carriers in the neutral back- 
ground (singly occupied sites) , can properly describe the 
Mott transition. 34 ~ 37 ' In particular in ref. 37, the D-H 
binding-unbinding mechanism is studied for a projected 
d-wave singlet state as well as a projected Fermi sea 
for the two-dimensional Hubbard model {t-t'-U model). 
A similar result was also reached using a Jastrow-type 
wave function. 38 ' Later, we have found in two dimensions 
that Mott transitions occur even in the wave functions 
in which a doublon must be necessarily accompanied by 
at least one holon in the nearest-neighbor (NN) sites. 
Namely, the state in which a doublon and a holon always 
tightly bind one another can be metallic (see §4.1). This 
finding requires a modification of the above picture 37 ' 
of Mott transitions through D-H binding and a simple 
release from it. 

In this paper, we address the following subjects by ap- 
plying a VMC method to the half-filled-band Hubbard 
model on the square lattice: (i) We extend the D-H at- 
tractive correlation factors so as to include some long- 
range types, and corroborate the decisive effect of D-H 
binding on nonmagnetic metal-insulator transitions, (ii) 
We introduce long-range D-D (and H-H) factors indepen- 
dently of the above D-H factors. Thereby, we can distin- 
guish the roles of the two factors for the Mott transition, 
(iii) We generalize the picture of the Mott transition so 
as to comprehend the one arising in the completely D- 
H bound state, and to treat it more quantitatively. By 
checking various wave functions, we are convinced that 
this conception is applicable to a wide range of systems, 
including Bose Hubbard models. 39 ' 

This paper is organized as follows: In §2, the method 
used in this paper is formulated. In §3, we discuss var- 
ious aspects of the VMC results. In §4, we propose an 
improved picture of Mott transitions, and confirm its ap- 
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plicability. Section 5 is assigned to summary. 

A part of the results in this study was published be- 
fore. 40 ) 

2. Formulation 

After brief introduction of the Hubbard model in §2.1, 
in §2.2, we describe the trial wave functions treated in 
this paper. In §2.3, we outline the VMC calculations. 

2.1 Hubbard model on square lattice 

We study the single-band Hubbard model with t,U ^ 
(3^24,41,42) wmcn j s fundamental to describe the physics 
of Mott transitions: 43 ^ 



< 



< 1) is a 



H = H t + Hu 



(1) 



where Ci„ is an electron annihilation operator of site i 
and spin a, di — n^riii and n ia = c\ a c ia . In this paper, 
we focus on the case of half filling on the square lattice 
with only the nearest neighbor (NN) hopping, 



H t = Ye k c\ a c k<y , 



6^ = —2t(cos k x + cos k y ) , 



(2) 
(3) 



and use t as the energy unit. In applying a variational 
Monte Carlo (VMC) method to this model, we use finite- 
size systems of L x L (= N s ) sites up to L = 18 
(7V S = 324) with the periodic and antiperiodic boundary 
conditions in x and y directions, respectively, to meet the 
closed-shell condition. 

2.2 Variational wave functions 

In §2.2.1, we explain the conventional NN D-H bind- 
ing wave function and the completely D-H bound state 
as its limiting case. In §2.2.2, we introduce a series of 
wave functions with long-range D-H binding factors and 
independent D-D (H-H) correlation factors. 

2.2.1 Nearest-neighbor correlation factor 

The simplest but fundamental trial wave function is 
the celebrated Gutzwiller wave function (GWF), 24 ) 



* G = Pg$f, 
where $f is the Fermi sea, and 

P G = J] [1 - (1 - g)dj] 



(4) 



(5) 



Here, g is the variational parameter which adjusts the 
doublon density, d = J2j(dj)/N s . It is known that GWF 

is metallic for U/t < oo. 28 ) To describe the Mott tran- 
sition, a wave function 31, 33 - 34 ) with a NN D-H binding 
correlation 44 ) has been often used: 



A(NN) 



nNN 



n (i hq : 



(6) 



where hj = (1 — nj^)(l — nji), \x (0 
variational parameter which controls the number of iso- 
lated doublons (D without H in its NN sites) and isolated 
holons, and f runs over the four NN sites of the site j. 
For fj, = 0, ^E'a(nn) is reduced to ^q- In the other limit, 
fi = 1, ^a(nn) becomes the completely D-H bound state, 



A(bind) 



Ptbind 
A 



$G=n{ i -^} cf>G - (8) 



In vE'A(bind): a doublon (holon) must be accompanied by 
at least one holon (doublon) in its NN sites, so that, su- 
perficially, ^Afbind) always seems insulating. However, it 
turns out to be metallic for small U /t, because a dou- 
blon (holon) can have multiple holons (doublons) in its 
NN sites, as will be discussed in §4.1. 

2.2.2 Long-range Jastrow factors 

According to the result of exact diagonalization in the 
one-dimensional half-filled-band Hubbard model, 33 ) the 
ground state of which is a paramagnetic insulator for any 
positive lift, 45 " 1 the magnitude of the coefficient of the 
basis having only one D-H pair with the D-H distance r 
decreases exponentially with r for large U ft. Assuming a 
similar situation arises in two dimensions, we introduce 
several types of long-range D-H attractive (A) correla- 
tion factors Pa, and the corresponding D-D and H-H 
repulsive (R) factors Pr: 



^a = II .Mlrfl) 



1- )[ (l-h j+ r) 



r£{rf} 

n (i 



dj+f) 



(9) 



Pr = 



n & (i*?o d. 



n u-^-ho 

fe{rf} 



+ hi 



1 n d 



ij+r) 



re{if} 



(10) 



where the index j of the outer product runs over all the 
sites, and ff (r* 1 ) indicates the vector from the site j 
to the nearest partner. Namely, we disregard more dis- 
tant ones, so that the index f of the products in eq. (9) 
runs over only the sites of the nearest D-to-H (in the 
first term) and H-to-D (in the second term) distances. 
This choice is reasonable in view of the strong-coupling 
expansion. 46 ) Correspondingly, we treat the product of 
f in cq. (10) in a similar way. To measure distances, we 
adopt the stepwise or Manhattan metric, in which, for 
instance, r — 2 for «-> (i + 1, j + 1), and the range 
is 1 < r < L. For each correlation, we consider the three 



Qi = d i IJ( 1 - h o+r) + ^ ]J(1 - d j+f ), (7) 
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forms: 



Mr) 
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1 



I — a exp 
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(a) exponential 

(b) power law 

(c) optimizing 



(11) 



1 



(a) exponential 

(b) power law 

(c) optimizing 



(12) 

where we fix /a(1) and /r(oo) at unity for (a) and (b). 



A(opt) 



Pl C) Pi C) *G, 



(15) 



'R(opt) - ' R 
'AR(opt) = -Pa 1 R 

where the superscripts of Pa and Pr correspond to the 
function types (a)-(c) in eqs. (11) and (12). For V^ar, we 
unify the function types of /a(^) and fn(r), for simplic- 
ity. In Table. I, we summarize the used wave functions. 

Note that the projector Pa Pr m e 1- (15) is different 
from the Jastrow factor used in related papers, 38, 47 ) es- 
pecially in that the magnitude of long-range D-D factor 
in them is connected to the inverse of corresponding D- 
H attractive factor, so that the D-D factor is necessarily 
repulsive as far as the D-H factor is attractive. 



(a) exponential type 




(b) power type 
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Fig. 1. (Color online) Weight of long-range D-H attractive corre- 
lations as a function of distance between D and H; (a) an expo- 
nentially decaying type [eq. 11(a)], and (b) a power-law decaying 
type [eq. 11(b)]. 



In eq. (11), each £ (£ r ) is a variational parameter. We 
consider two typical forms: (a) exponentially decaying 
and (b) power-law decaying, as shown in Fig. 1. Fur- 
thermore, in (c), we do not assume a specific form of 
/a(?")> and optimize all £ r 's (2 < r < L) simultaneously 
as variational parameters. The type (c) is the best form 
of f\{r), here; the other types are specific cases of (c). 
In eq. (12), we again assume the repulsive correlation 
becomes weaker, as r increases. Corresponding to /a(0, 
we consider three forms for fB,(r): (a) an exponentially 
and (b) a power-law decaying types, both with the pa- 
rameters a adjusting the weight of r = 1, and /3 con- 
trolling the decaying length. In (c), we optimize all a r 's 
(1 < r < L) simultaneously. Notice that here we permit 
a r to be larger than 1, meaning Pr possibly works as 
D-D (H-H) attractive correlations. 

In this work, we study a series of wave functions by 
combining the above correlation factors as follows: 



A(cxp) 
R(cxp) 



AR(cxp) 



Pf^G, 

Pr^G, 
p(a) p(a)^rv 



A(pow) 
R(pow) 



p{ b) * 
p£>* 



AR(pc 



(13) 



(14) 



Pl b) Pi b) *G, 



2.3 Variational Monte Carlo calculations 

In this subsection, we briefly describe the outline of 
VMC calculations implemented in this study. 

A correlated measurement or optimization- VMC tech- 
nique 48 ) is used to optimize variational parameters up 
to 2L. In the non-linear minimization process of energy 
expectation values for the wave functions with many pa- 
rameters, we adopt a quasi-Newton method, in which 
gradient vectors are effectively calculated by recently 
proposed formulae, 49 ^ and Hessian matrices are approxi- 
mated by Broyden-Flecher-Goldfarb-Shanno formula, 50 ^ 
the use of which does not affect the exactness of optimiza- 
tion itself. In coding, we refer to an algorithm offered by 
Ibaraki and Fukushima. 51 -' For wave functions with a few 
parameters, we use a simple linear optimization together. 

In both algorithms, parameters as well as energy con- 
verge typically after first several rounds of iteration with 
different fixed sample sets; in each set we generate typ- 
ically 2.5 x 10 5 particle configurations using Metropo- 
lis algorithm. After the convergence, we continue excess 
rounds (10-20 times) of iteration in the optimization pro- 
cess with successively renewed configuration sets. We de- 
termine the optimized values by averaging the data ob- 
tained in the excess rounds; in averaging, we exclude 
scattered data beyond the range of twice the standard 
deviation. The optimal value is an average of substan- 
tially more than several million samples. The variational 



Table I. Summary of trial wave functions. In the second column, 
we abbreviate the type of correlation, with GW being the onsite 
(Gutzwiller) repulsion. 





correlation 


correlation 


param. 


abbreviation 


type 


range 


number 


GWF 


GW 


onsite 


1 


A(NN) 


GW+DH 


NN 


2 


A(bind) 


GW+DH 


NN 


1 


A(cxp) 


GW+DH 


long 


2 


R(exp) 


GW+DD 


long 


3 


AR(exp) 


GW+DH+DD 


long 


4 


A(pow) 


GW+DH 


long 


2 


R(pow) 


GW+DD 


long 


3 


AR(pow) 


GW+DH+DD 


long 


4 


A(opt) 


GW+DH 


long 


L 


R(opt) 


GW+DD 


long 


L+l 


AR(opt) 


GW+DH+DD 


long 


2L 
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energy and significant parameters [g, /i, £, etc.] are ob- 
tained with sufficient accuracy, but the determination of 
insignificant parameters [/a(»") and /r(V) with r > 8, 
etc.] is difficult, because E depends on them only very 
slightly, in other words, particle configurations determin- 
ing them appear extremely rarely. Anyway, such parame- 
ters have little influence on E and other quantities. Phys- 
ical quantities are calculated typically with 2.5 x 10 5 re- 
newed configurations generated by the optimized param- 
eter sets. 

Since in Mott critical regimes, the global minimum 
becomes more competitive with other local minima as L 
increases, accurate energy minimization sometimes be- 
comes not easy. This leads to the scattered data points 
near U c /t in some figures. 

3. Results 

In this section, we discuss the results of VMC calcu- 
lations. In §3.1 and §3.2, energies and other quantities 
obtained by various wave functions are studied, respec- 
tively, in view of the Mott transition. In §3.3, the system- 
size dependence of the Mott critical values are discussed. 
In §3.4, the differences among various D-H factors are 
studied. In §3.5, we discuss the effects of repulsive long- 
range factors. 

3. 1 Behavior of energy 




4 8 12 jy /f 16 20 24 28 



Fig. 2. (Color online) Comparison of variational energies among 
various trial wave functions as function of correlation strength for 
L = 16. The pale solid line is a guide line proportional to —t/U 
expected from the strong-coupling expansion. The inset shows a 
magnification near the crossing point of GWF and A (bind). 



First of all, let us discuss the behavior of total energy 
per site E/t. In the main panel of Fig. 2, we compare 
E/t among various wave functions in a wide range of 
U/t. Here, we notice GWF and the completely bound 
state ^E'A(bind)- It is known that GWF is metallic for 
U/t < oo and a relatively good for U sufficiently smaller 
than the critical value of the Brinkman-Rice transition, 
Ubr = 12.97i. 25 ' On the other hand, \&A(bind) is insu- 
lating for U/t ^ 3 (see §4.1), which fact is supported 



by the coincidence with the result of strong coupling ex- 
pansion oc —t/U, 4 ^ as shown in Fig. 2. E/t of GWF is 
much lower than that of ^A(bind) f° r U/t < 7.4, while the 
relation is reversed for U/t > 7.4, suggesting the phase 
switches from metal to insulator at U/t ~ 7.4. Now, we 
look at the D-H binding wave functions (^a)- E/t of all 
^a's [-E'(A)] except for V&Afbind) behave similarly, namely, 
E(A) is very close to £(GWF) for U/t Sa 5, somewhat 
smaller than both £(GWF) and -E(^A(bind)) f° r the in- 
termediate values of U/t, and approaches ^(^Afbind)) 
for U/t ^ 9. Hence, a wide class of D-H binding wave 
functions probably induces Mott transitions at U ~ W 
\W(— 8t): band width]. This is consistent with the previ- 
ous result for a short-range D-H binding wave function, 
in which a first-order Mott transition occur at C7 c /t=8.59 
and 8.73 for L = 16 and 18, respectively. 37 ^ 

The detailed behavior of E(A) for intermediate U/t 
is different among different types of D-H factors, as 
shown in the inset of Fig. 2. For example, -E^Atexp)) 
exhibits a clear cusp at U ~ 8.12£ (= U c ), and has a 
higher value than those of other E(A)'s for U < U c , 
whereas i?( x I , A(opt)) exhibits smooth behavior in this 
regime, and has a lower value. The energy of the best 
function x I , AR(opt) is broadly similar to -E/(^A(opt)) f° r 
U < U c , but has an appreciably lower energy for U > U c . 
We will return to this subject in §3.4. 




Fig. 3. (Color online) The expectation values of (a) kinetic en- 
ergy and (b) doublon density (interaction energy) are compared 
among the identical wave functions treated in Fig. 2. The insets 
show magnifications near the Mott critical points. 
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Next, we consider the components of energy. Although 
the total energy should be a continuous function of U /t, 
its components can exhibit stronger critical anomalies at 
U c /t if the transition is first order. Figures 3(a) and 3(b) 
shows the kinetic energy, E^ in = (H t )/N s , and doublon 
density (substantial interaction energy), d = E- lnt /U — 
(Hu)/(UN S ), respectively. As evidently seen in the in- 
sets, 5 , A(cxp), which exhibits a clear cusp in E/t, has 
discontinuities at U c /t in E^in and d. This is an obvious 
sign of a first-order transition; actually we have observed 
a hysteresis around the critical point. Other wave func- 
tions behave more mildly at this system size, but the be- 
havior evolves into more first-order-like as L increases, as 
we will discuss in §3.3. The fact that Ek- ln (d) of every D- 
H binding wave function abruptly increases (decreases) 
in the critical region as U/t increases is consistent with a 
requisite of the Mott transition: The metal-to-insulator 
transition is driven by reducing the interaction energy at 
the cost of the kinetic energy. This is in sharp contrast 
with antiferromagnetic and superconducting transitions 
arising in strongly correlated regimes. 21,35 ) 

3.2 Critical behavior of physical quantities 

To corroborate the realization of Mott transition in 
and estimate the critical value more accurately, we 
consider other physical quantities. 




xo,o) £ M(x,7t) r(o,o) 

Fig. 4. (Color online) The momentum distribution function is 
shown along the path, (0, 0)— »(-7r, 0)— >(tt, it)— >(0, 0), calculated 
with the wave funcition AR(opt) for various values of U/t near 
Uc/t ~ 8.27 for L = 16. The vertical dash-dotted lines indicate 
the positions of fep in the metallic cases. 

First, we discuss the momentum distribution function, 

(16) 

In Fig. 4, we show, as a typical case, the result of 
\E f AR(opt) ; which exhibits first-order-like critical behavior 
in energy at U c /t ~ 8.30 (Fig. 3). Let us pay attention to 
the behavior around the Fermi surface near k = (7r,0). 
For small values of U/t (< 8.3), a discontinuity of n(k) 
is evident at fcp, whereas the magnitude of discontinu- 
ity abruptly becomes small at U ~ U c and remains very 



small for large U/t (> 8.3). To discuss quantitatively, we 
actually measure the jump of n(k), namely quasiparticle 
renormalization factor Z, at the X point (jr, 0): 

Z = n~(k)\ k ^ X -o - n + (k)\ k ^ x +o, (17) 

where n~(k) [n + {k)\ is the fitting function of the seg- 
ment F-X [X-M] of n(k) given by the third-order of 
least squares method. According to the Fermi liquid the- 
ory, Z becomes zero for insulating states. As shown in 
Fig. 5, Z of every D-H binding function almost vanishes 
at U/t — 8.2-9.0, although the detailed behavior some- 
what differs among different wave functions. The small 
residual values for large U ft are owing to the finite sys- 
tem sizes. Thus, we may regard the states for U/t ^ 9 
as insulating. As for vE f AR(opt)i % most markedly drops at 
8.2-8.3, which coincides with U c /t evaluated from energy. 
In contrast, Z of a metallic state GWF asymptotically 
approaches zero, as known. 34 ) 

Second, we consider the charge density correlation 
function in the wave-number space, 

N @ = w a ? eWi ^° m - n2 > ( 18 ) 

which is known, within the variation theory, to behave 
as N(q) oc \q\ for \q\ — > unless an excitation gap opens 
in the charge degree of freedom, whereas N(q) oc \q\ 2 if 
a charge gap opens. Figure 6(a) shows N(q) of ^ARfapt)! 
the behavior for small values of \q\ abruptly changes be- 
tween U/t = 8.25 and 8.3, which again coincide with U c /t 
determined by other quantities. For U > U c , the behav- 
ior becomes |<7| 2 -like, suggesting a charge gap opens. 
Third, we study the spin correlation function, 

S(q) = ^- Z Z emfi ^ ) ( S i S ])- (I 9 ) 

S ij 

In Fig. 6(b), S(q) of "J/AR^pt) is plotted for the same 
range of U/t as in (a). For small the behavior is al- 
ways S(q) oc \q\, suggesting the low-lying spin excita- 




U/t 



Fig. 5. (Color online) Quasiparticle renormalization factor as 
function of U/t, estimated from discontinuities in n(k) at k = 
(■7T, 0) for various type of wave functions. The inset shows the 
magnification near the Mott critical points. A(bind) will be dis- 
cussed later. 
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tion is gapless both in metallic and insulating phases. 
On the other hand, the magnitude at the AF nesting 
vector q = (tt,w) increases as U/t increases, in particu- 
lar, markedly at U/t = 8.3. Thus, the D-H binding wave 
function shows strong inclination toward the antiferro- 
magnetic order, especially in the insulating case. This 
aspect is the same as the short-range D-H binding wave 
functions. 37 ' 

3.3 System-size dependence of critical point 

In this subsection, we consider the system-size depen- 
dence of the Mott critical point, which we have disre- 
garded to this point. 

The optimized variational parameters controlling the 
D-H binding strength is important to definitely deter- 
mine the critical values. In the main panel of Fig 7(a) , the 
optimized 1 /£ in the wave functions of the exponentially- 
decaying type [eq. (11a)] is plotted for three system 
sizes. For small systems like L = 10, the variation of 
l/£ is smooth and clear critical behavior is not seen, 
whereas systems with L > 14 exhibit clear discontinu- 
ities. Such tendency is often observed when finite sys- 
tems are used, 52 ) because the phase transition is well 
defined for L = oo. As another example, we show the 
result for the short-range D-H correlation in the inset of 
Fig. 7(a). 37 - 39 ) 

It is also a general tendency of the D-H binding wave 
functions that the Mott critical value increases as L in- 
creases. 37,39 ' This stems from the great system-size de- 
pendence in E/t in the insulating side of the transition 
point, compared with in the metallic side. This is re- 
flected in some decrease of the critical value by adding 
repulsive Jastrow factors [ x t , AR(oxp)]i which improve E/t 
especially in the insulating side, as in the inset of Fig. 2. 




(0,0) (a;0) q {k,7z) (0,0) 

Fig. 6. (Color online) (a) Charge density structure factor of 
AR(opt) for some values of U/t near U c /t ~ 8.3. (b) Spin struc- 
ture factor of AR(opt) for the same condition as (a). 



These tendencies are common to other D-H binding 
wave functions; as another example, in Fig. 7(b), we show 
optimized £ in the power-law decaying correlation factor 
cq. (lib). In x I , A(pow) J the critical behavior is less clear. In 
case the Mott critical point cannot be clearly determined, 
we estimate U c /t from the maximum of the decreasing 
rate of doublon density d, 

dd = d[(U + AU)/t] - d[U/t] 
d(U/t) AU/t ■ 1 ' 




U/t 

Fig. 7. (Color online) The behavior of optimized variational pa- 
rameters which controls D-H binding strength is compared 
among three system sizes near the Mott critical points. The re- 
sults for three types of wave functions are shown: (a) l/£ for 
A(exp) and AR(exp). The inset shows fi for A(NN). (b) £ for 
A(pow) and AR(pow). 



Table II. The Mott critical values, Uc/t, estimated from eq. (20) 
is entered in the first line of each wave function. The results of 
five system sizes are compared. The figures in the second line 
represent the critical values determined by the crossing point of 
{dh and £dd, which will be discussed in §4. The figures with 
* for A(bind) are estimated from Z, which will be explained in 
§4.1. 





L = 10 


L = 12 


L = 14 


L = 16 


L = 18 


A(NN) 


7.85 


8.25 


8.475 


8.575 


8.675 


'dh = ^DD 


8.85 


8.725 


8.575 


8.575 


8.875 


A(exp) 


7.75 


7.875 


8.025 


8.125 


8.275 


*DH = ^DD 


7.95 


7.875 


8.025 


8.125 


8.275 


AR(cxp) 


7.65 


7.775 


7.925 


8.025 


8.125 


*DH = ^DD 


7.75 


7.775 


7.925 


8.025 


8.125 


A(pow) 


7.55 


7.95 


8.375 


8.575 


8.675 


*DH = ^DD 


8.55 


8.475 


8.525 


8.625 


8.725 


AR(pow) 


7.65 


7.95 


8.175 


8.375 


8.525 


'DH = &D 


8.35 


8.275 


8.325 


8.425 


8.525 


A (opt) 


7.15 


8.075 


8.275 


8.525 


8.875 


*DH = ^DD 


8.75 


8.775 


8.825 


8.925 


9.025 


AR(opt) 


7.45 


7.925 


8.125 


8.275 


8.375 


*DH = ^DD 


8.225 


8.125 


8.175 


8.275 


8.375 


A(bind) 


2.10* 


2.28* 


2.365* 


2.39* 


2.42* 


^DD = ^ 


2.00 


2.175 


2.275 


2.375 


2.425 
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The Mott critical values thus obtained are summarized 
in Table II. 

3.4 Effects of different D-H binding factors 




Fig. 8. (Color online) The appearance rates of nearest (a) D-H 
and (b) D-D distances of r are plotted for r = 1-7 as a function 
of U/t. Both Wdh and Wqd are calculated with AR(opt). The 
vertical dashed lines indicate the Mott critical point. 

First of all, we look at the distribution of distance from 
a doublon (holon) to the nearest holon (doublon), which 
is denoted simply by r, here. We indicate the appearance 
rate of r by Wdh(p)) which satisfies, 

L 

^W DH (r) = l. (21) 

r=l 

In Fig. 8(a), Wdh(»") with r < 7 for the best function 
^AR(opt) is depicted versus U/t. In the metallic regime, 
Wdh(1) is predominant, but Wdh(^) with r > 2 has an 
appreciable weight especially at U/t ~ 6; there, a dou- 
blon is detached from holons to some extent. Meanwhile, 
in the insulating regime, the weight is almost concen- 
trated on r = 1, indicating D-H pairs are confined within 
mutually NN sites. In a similar way, we define WddM as 
the appearance rate of doublons (holons) with the D-to- 
D (H-to-H) distance of r. As shown in Fig. 8(b), Wnn(r) 
of large r increases and Wdd(1) and Wdd(2) rapidly 
decreases, as U/t increases. In the insulating phase, dou- 
blons tend to keep away from each other. 

In the following, we discuss the effects of different D- 
H binding factors P A X ' (x = a, b, c, NN), without in- 
troducing repulsive correlations. Figure 9 compares the 
optimized weights fh{r) [eq. (11)] among the four D-H 
projectors. The plotted data are restricted to what sat- 
isfies the following condition: 

p{r) > Pmin, (22) 

where r has the same meaning as in the preceding para- 
graph, and p(r) [= WdhO") x d ] denotes the appearance 
probability of the doublon (holon) of r. After a search, 
we put p m j„ = 4 x 10~ 4 , which corresponds to the prob- 
ability that doublons or holons of certain r appear 10 4 
times in 2.5 x 10 5 samples for the system of L = 10. If the 
condition (22) is not satisfied for r = r* , the optimized 



/a(?-*) become statistically unreliable in the present cal- 
culations, and corresponding samples actually make only 
an imperceptible contribution to the averages. 

We start with the analysis of the optimizing-type wave 
function ^/A(opt)- F° r any values of U/t, fx{r) rapidly 
decreases for r Si 3, but becomes almost constant for 
r *1 3. Note that this long-range behavior of /a(^) is 
convenient for the conductive nature in the metallic 
regime (U < U c ), namely, a doublon is released from 
the bondage of holons, once a doublon goes three lattice 
constants away from holons. The effective range of /a(^) 
satisfying the condition (22) becomes the widest (r 8) 
near U c /t, but abruptly shrinks in the insulating phase, 
as expected from the data in Fig. 8(a). In addition, for 
U > U c , the magnitude of /a(^) itself is small. As a re- 
sult, a doublon and a holon come to confine each other in 
a narrow D-H pair domain. We will return to this topic 
as to the mechanism of Mott transitions in §4. The pro- 
priety of other 'J a depends on how properly /a(^) of "J a 
can imitate that of v I , A(opt) ■ 

Regarding a short-range D-H factor, the reason why 
simple \E'a(nn) is unexpectedly good for U < U c (inset 
of Fig. 2) is that /a(0 of ^a(nn) is constant for r > 2 
and resembles that of ^Afopt) except for /a(2), as shown 
in Fig. 9 for the data of U/t = 7. On the other hand, 
the reason why E/t of ^a(nn) is relatively high in the 
insulating regime, as compared with the other ^a, stems 
from an underestimate of /a (2), as seen for U/t = 12. 

We turn to the long-range D-H factors. In the metallic 
regime, /aW s of * A ( P ow) and especially of * A (exp) are 
underestimated for small r and overestimated for large 
r, to be optimized as a whole. Consequently, ^/A(exp) has 
appreciably higher energy than other V^a's; this behavior 
lowers the Mott critical value of 'I'A(exp)- On the other 
hand, in the insulating regime, /a(»*)'s of both "^(cxp) 
and v I / A(pow) almost coincide with that of ^A^pt^ be- 
cause the long-range part of /a(»") (r ^ 3) substantially 
vanishes. As a result, -E( v I , A(oxp)) and -E( x I , A(pow)) be- 
come as good as i?( x I , A(opt))- 




Fig. 9. (Color online) Comparison of optimized D-H attractive 
correlation weight /aM among four D-H binding wave functions, 
A(NN), A(exp), A(pow) and A(opt), for various values of U/t. 
We omit the data which do not meet the condition (22). 
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3. 5 Effects of repulsive intersite correlations 

First of all, we compare contributions to the improve- 
ment of energy between D-H attractive and D-D (and 
H-H) repulsive correlation factors. In Table III, the op- 
timized total energies of GWF and the optimizing-type 
wave functions eq. (15) are summarized. The results of 
other types of wave functions eqs. (13) and (14) are ba- 
sically identical as far as the effect of repulsive factors is 
concerned. The wave function 'I'R(opt) [ ec l- (15)], in which 
only the repulsive factor Pr [eq. (10) with eq. (12c)] 
is applied to GWF, improves E/t only very slightly on 
GWF for any value of U/t. Furthermore, ^'R(opt) never 
induces a Mott transition, like GWF. Thus, the repul- 
sive correlation Pr, by itself, make no substantial im- 
provement on GWF. In contrast, as already discussed 
for Table II and Fig. 2, v I / A(opt) induces a Mott transi- 
tion, and make a great improvement in E/t on GWF 
for U ~& U c , conclusively showing that the essence of 
Mott transition is included in the D-H binding corre- 
lation. * A R(opt) (= Pi*A(o P t)) further reduces E/t for 
U > U c - Although the decrement in energy thereby is as 
small as 2% of P(*A( op t)) [U/t = 9 and 12], the addition 
of Pr shifts the Mott critical point to a somewhat lower 
value (Table II), and make the critical behavior more 
first-ordcr-like (inset of Fig. 5). 



Table III. Comparison of variational energies E/t among GWF 
and the optimizing type of wave functions, eq. (15), for L = 16. 
The digits in the brackets denotes the error in the last digits. 



U/t 


GWF 


R(opt) 


A(opt) 


AR(opt) 


1.0 


-1.3865(2) 


-1.3866(3) 


-1.3867(4) 


-1.3867(3) 


7.0 


-0.374(1) 


-0.375(2) 


-0.418(2) 


-0.418(2) 


7.5 


-0.323(2) 


-0.323(2) 


-0.379(2) 


-0.379(2) 


8.0 


-0.276(2) 


-0.277(2) 


-0.348(2) 


-0.348(3) 


8.5 


-0.234(2) 


-0.234(2) 


-0.325(3) 


-0.330(2) 


9.0 


-0.196(2) 


-0.197(2) 


-0.310(2) 


-0.317(2) 


12.0 


-0.061(2) 


-0.063(2) 


-0.255(3) 


-0.261(2) 



0.02 




12 , 16 

U/t 



Fig. 10. (Color online) The differences of total, kinetic and in- 
teraction energies between AR(opt) and A(opt) are plotted as a 
function of U/t for L = 16. If the value is negative, AR(opt) has 
a lower energy than A(R). The large deviations near U c /t should 
be neglected, which stem from the discordance of U c /t. 



Next, we consider the differences of kinetic and of in- 
teraction energies between Vl'AR(opt) and ^A(opt) : 

AP r = P r [AR(opt)] - P r [A(opt)] (23) 

where the suffix T denotes "kin", "int" or "tot". In 
Fig .10, the three kinds of AP are plotted versus U/t. 
In the metallic regime, not only AP to t remains zero, but 
also both APkin and APi nt are zero except for irreg- 
ular accidental deviations. Thus, Pr does not modify 
■E(*A(opt)) for U < U c . Nevertheless, for U/t > 10 in 
the insulating regime, AP^n exhibits appreciable neg- 
ative values, according to AP tot . Inversely, AP pot has 
regular positive values; Pr with the aid of Pa reduces 
Pkin at the cost of E- lnt in the insulating phase. Thus, Pr 
compensates for the excess of D-H binding effects. 



1.2 



1.1 



^1 



0.9 



1 


Y 


(a)l-g " 


1= 16 AR(opt) 




<b)Mr) 


- 2 - 




- — U c /t=S.21- 



0.9 



U/t 



Fig. 11. (Color online) (a) Optimized onsite repulsive (Gutzwil- 
ler) parameter g near the Mott critical point, (b) Optimized long- 
range repulsive correlation weight fn(r), eq. (12c). The results in 
(a) and (b) are obtained in an identical calculation for AR(exp). 



Finally, we discuss the repulsive correlation weight 
for ^AR(opt)j m which /r(7") is optimized for each r. 
In Fig. 11(b), fn(r) for r < 4 is shown near U c /t. 
In the metallic regime, fn(r) is almost unity, although 
there is some fluctuations, indicating the repulsive inter- 
site correlation is virtually ineffective. Correspondingly, 
^AR(opt) seldom improves E/t on that of ^A(opt) f° r 
U < U c (Table III). Meanwhile, in the insulating regime, 
only /r(2) is slightly lowered from unity, causing the 
sudden drop of W^ D d(2) for U > U c [Fig. 8(b)]. 53 ) The 
behavior of fn(r) > 1 for r > 3 indicates that D-D and 
H-H correlations are attractive rather than repulsive in 
a long-range part. 54 ) This behavior promotes the repul- 
sion between doublons in proximity, as seen in Fig. 8(b), 
where Wud(t) for r > 3 is larger than Wdd(1)- These 
effects of fn(r) causes a small but steady improvement in 
energy on ^A(opt) (Table III). Incidentally, the scattered 
data, especially for U < U c , stem from the complemen- 
tarity between the onsite and intersite repulsive correla- 
tions. The behavior of g in Fig. 11(a) is quite opposite 
to that of /r(t") in Fig. 11(b). Owing to the balance be- 
tween the two, the resultant physical quantities become 
smooth, for instance, d in Fig. 16. In this point, the pa- 
rameter space has redundancy. 
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4. Improved Picture of Mott Transitions 

In §4.1, we discuss the Mott transition arising in the 
completely D-H bound state, and provide a reformed pic- 
ture of Mott transitions to comprehend it. In §4.2, we 
show that this picture is applicable to various D-H bind- 
ing wave functions. 

4--1 Mott transition in completely D-H bound state 

In §3.1, we supposed that the completely D-H bound 
wave function ^Afbind) [ ec i- (8)] represents a typical in- 
sulating state, but this supposition is not true for small 
values of U/t. The total energy of 'il'A(bind) m Fig. 2 ex- 
hibits an abrupt change of curvature at U/t ~ 2.5. Cor- 
roboratively, as shown in Fig. 5, the quasiparticle renor- 
malization factor Z has finite values for U/t Ss 2.5. Thus, 
it is certain that ^A(bind) exhibits a Mott transition and 
is metallic for small U/t. Then, we estimate the Mott 



critical values of 



A(bind) 



from the extrapolation of Z 



for U/t 2 using the third-order least squares method, 
as shown in Fig. 12, and summarize them in Table II. In 
*A(bind)! doublons (holons) must be accompanied by at 
least one holon (doublon) in its NN sites. The metallic 
state satisfying this condition cannot be represented sim- 
ply by the unbinding of D-H pairs proposed in ref. 37. In 
what follows, we consider a microscopic picture of Mott 
transitions which comprehends the case of ^AOind) ■ 




U/t 



Fig. 12. (Color online) The Mott critical value U c /t for com- 
pletely D-H bound state A(bind) is estimated in two ways. The 
vertical dashed lines indicate the values estimated from the ex- 
trapolation (Z*) to zero of the quasiparticle renormalization fac- 
tor Z. The other estimate is given by the crossing point of £p H 
(= 3) and ££ D (= 1/V~d). The results are listed in Table. II. 



For simplicity, we take up a one-dimensional case, 55 - 1 
and first assume that U/t is small and the doublon (and 
holon) density is sufficiently high. Actually, for U/t ^ 1.5, 
d of v I / A(bmd) becomes higher than those of the other D-H 
binding states [see Fig. 3(b)]. Then, a doublon are often 
accompanied by multiple holons in the NN sites, so that 
it can propagate independently of holons as a plus charge 
carrier, satisfying the complete binding condition: 

[• • DHD^HDH •■]->[•• DH^HD^DH • •] -> 

[■ • D^HHDD^H ■■]->[■ ■<-> HDHDHD^ • ■] -> (24) 



Here, pay attention to the dotted D and H, which propa- 
gate by exchanging the positions according to the double- 
headed arrows. Thus, the state becomes metallic. On the 
other hand, when U/t becomes large and the doublon 
density becomes sufficiently low, most D-H pairs become 
mutually detached, as [• • f! HD |f HD || • •]• Then, a 
doublon becomes unable to propagate independently of 
holons, and vice versa, satisfying the completely bound 
condition. Consequently, charge fluctuation is confined 
locally, resulting in an insulator. Summing up, the con- 
duction in a metallic state requires that D-H pairs should 
contact with one another in sequence; the Mott critical 
value can be specified by U/t at which D-H pairs are 
mutually detached. 

From the above argument, we find it convenient for 
general discussions of Mott transitions to introduce two 
characteristic length scales, the D-H binding length £qh 
and the D-D exclusion length Idd , which are generally a 
function of U /t. We postulate that the attractive corre- 
lation factor Pa produces D-H pairs of a binding length 
£dh according to U/t; £dh approximately corresponds to 
the size of a D-H-pair domain, in which at least one dou- 
blon and one holon must exist. £dd broadly represents 
the distance between two D-H pairs. Using Idh and €dd, 
we give a microscopic picture of Mott transitions, which 
is schematically shown in Figs. 15(c) and (d) for general 
D-H binding wave functions. In the insulating phase, the 
relation ^dh < ^dd holds, indicating that the domains 
of D-H pairs do not usually overlap, at least, not in se- 
quence. Consequently, most D-H pairs are isolated and a 
doublon and a holon are confined within ^dh, resulting in 
only local charge fluctuation. To this point, the picture is 
basically identical with the previous one. 37 - 1 In the con- 
ductive phase (U < U c ), £bh becomes larger than ^dd, 
indicating the domains of D-H pairs overlap with one 
another. Then, a doublon in a D-H pair can exchange 
a partner holon with a holon in an adjacent D-H pair, 
when the two holons are in the overlapped area. As a re- 
sult, a doublon and a holon can move independently as 
charge carriers by successively exchanging the partner. 
As U/t is varied, a Mott transition takes place when £dh 
becomes equivalent to ^dd, which is roughly 1/y/d and 
generally a monotonically increasing function of U/t. In 
this framework, it is of primarily important to determine 
^dh and ^dd appropriately. 

Now, we apply the above framework to the special case 
of ^A(bind)- We postulate that a domain of a D-H pair 
for l &A(bind) consists of 3 x 3 lattice sites, as shown in 
Fig. 13. Note that the size of this domain for \I/A(bind) 
is constant, irrespective of U/t, because a doublon and a 
holon (or holons) are tightly bound in the nearest neigh- 
bor site(s). For ^dh < 3, some domains of D-H pairs 
mutually overlap, and sequence like DHDHDH in (24) 
possibly appears, whereas for foH > 3 the domains of 
D-H pairs separate from each other on average. Thus, it 
seems appropriate to put £p H = 3. Here, we add a star to 
distinguish it from the form for ordinary (not completely 
bound) ^a's discussed later. As for the D-D exclusion 
length, we simply put ^ D = 1/ yd, as mentioned above. 
In Fig. 12, we plot £^ D versus U/t. The values of U/t 
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at which 



DH 



= £ DD are summarized in Table. 



II, and 



approximately consistent with the Mott critical values 
estimated from Z . Thus, the above picture with a broad 
estimate of ^dh and £dd yields a justifiable result to this 
special case. 

4-2 Application to ordinary D-H wave functions 

Here, we study whether the above picture of Mott tran- 
sitions is applicable to ordinary D-H binding (not rigidly 
bound) wave functions V&a and to those with repulsive 
factors ^ar- In contrast to Vl/A(bind)i the distance from 
a doublon to its nearest holon(s), (t*dh)) varies in the or- 
dinary V^a and \&ar, so that not only Idd but also Inn 
should depend on U/t. In Fig. 14(a), we show the be- 
havior of the averages of the nearest D-to-H and H-to-D 
distances (run) and of the nearest D-to-D and H-to-H 
distances (ruu) calculated with ^t'AR(opt) are plotted as 
a function of U ft. Figure 14(b) shows the standard devi- 
ations of (r A ) (A = DH or DD), 



OA 



M 



(25) 



£ = 3 

1 DD J 





! 


! ! 




ci'ix : 


: xt'd 






I 



Fig. 13. (Color online) Two domains of D-H pairs for the com- 
pletely bound state A(bind) (domain size £dh is 3) are schemat- 
ically shown at the Mott critical case: £fjd = ^DH- Doublons 
and holons are represented by circles with and without arrows, 
respectively. The empty squares indicate singly occupied sites. 




U/t 



5 u/t 10 



Fig. 14. (Color online) (a) The nearest D-to-H and H-to-D dis- 
tances and the nearest D-to-D and H-to-H distances calculated 
with AR(opt) are plotted for three system sizes as a function of 
the interaction strength, (b) Standard deviations of (rrjij) an <3 
of (tub) shown in (a). 



where the index i runs over all the doublons and holons 
in all the measured samples, and M indicates their total 
number. a\ behaves similarly to (r" A ). Allowing for the 
meaning of ^dh and £dh, we put, generally, 



«DH = V DH/ + CTdH; 
£dD = (n3D> - ODD- 



(26) 
(27) 



Namely, £dh broadly represents the maximum radius of a 
D-H pair domain, over which a doublon does not separate 
from holons, as in Fig. 15(a), and ^dd the minimum D- 
D length, under which doublons do not approach each 
other, as in Fig. 15(b). 




' Mott trans. \ 
(' dh > I dd (metal) £ DH a £ DD £ DH < t DD (insulator) 

Fig. 15. (Color online) Schematic figures of microscopic mecha- 
nism of the Mott transition, (a) The solid circle of radius £fjh 
denotes the domain in which a holon can itinerate, when the 
partner doublon is located at the center, (b) The solid circle of 
radius £dd denotes the forbidden area where a doublon cannot 
enter, when another doublon is situated at the center. In (c) and 
(d), the behavior of doublons (D) and holons (H) is illustrated 
for metallic and insulating states, respectively, in ordinary D- 
H binding wave functions. The explanation is given in §4.1. The 
two phases are distinguished by comparing the magnitude of ^dh 
and <?dd- 



In the metallic regime, (run) and ctdh gradually in- 
crease as U/t increases as shown in Fig. 14, chiefly be- 
cause the densities of doublon and holon are reduced by 
U/t [see Fig. 3(b)]. This effect exceeds the D-H binding 
effect of Pa- Consequently, a doublon somewhat separate 
from holons. At the Mott critical point, however, (tdh) 
and ctdh suddenly drop, and asymptotically approach 1 
and 0, respectively, in the insulating regime, owing to 
the predominant D-H binding effect. In contrast, (rno) 
and (Tdd monotonically increase as U /t increases, owing 
to the steady decrease of d. In Fig. 16, £dh and Idd ob- 
tained through eqs. (26) and (27) are plotted. The Mott 
critical points U c /t are estimated from the steepest de- 
scent of an order parameter of the Mott transition d, 
and are listed in Table II. It is found that £dh intersects 
^dh almost at U c /t. Namely, the Mott critical point is 
estimated also from the condition, 



-DH 



('DD* 



(28) 
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The values of U c /t thus obtained are listed in Table II. 
Thus, the picture of the Mott transition introduced in 
§4.1 [Figs. 15(c) and 15(d)] is applicable to x I , AR(opt)- 




U/t 



Fig. 16. (Color online) The D-H binding length <dh and the D- 
D exclusion distance £dd for AR(opt) obtained from the data 
in Fig. 14 are plotted as a function of U/t. For comparison, the 
doublon density and the critical values obtained thereby (vertical 
dashed lines) are added for the same systems. 

Now, we confirm whether the above scheme with 
eqs. (26) and (27) are effective also for other wave func- 
tions. We have made similar analyses for various types 
of \&a and ^ar in Table I. As a results, the Mott critical 
values determined under the condition eq. (28) coincides 
with U/t estimated from d for every wave function with 
sufficient accuracy at least for large L. As an example, 
in Fig. 17, we plot the same quantities as in Fig. 16 for 
four representatives. 
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Fig. 18. (Color online) D-H binding length £dh and D-D exclu- 
sion distance £fjd calculated with Gutzwiller wave function as 
function of interaction strength. 

Finally, we touch on the wave functions with only re- 
pulsive correlation factors, namely GWF and the three 
\I/r's in Table I. As mentioned in §3.5, these wave func- 
tions do not undergo Mott transitions, and are always 
metallic for U/t < oo. In Fig. 18, we plot ^dh and £dd of 
GWF versus U/t, as a typical example, £dd is a mono- 
tonically increasing function of U jt in the similar way as 
"J A and "Jar, whereas €dh is also a monotonically in- 
creasing function and is always about twice larger than 



£dd- Namely, the relation £dh > &d always holds and 
^dh never crosses £dd- The situation is the same for the 
three 'I'r's. Thus, the new picture on the Mott transition 
properly works for the repulsive Jastrow-type wave func- 
tions. From this argument, it is certified again that the 
D-H binding effect is the essence of the Mott transition. 

5. Summary 

In this paper, we have studied the nonmagnetic Mott 
transition in the Hubbard model on the square lattice, 
using a variational Monte Carlo method. In the trial wave 
functions, we introduce long-range doublon-holon attrac- 
tive and doublon-doublon (holon-holon) repulsive corre- 
lation factors, in addition to the onsite repulsive factor. 
We recapitulate the main results below. 

(1) We confirmed that the D-H binding correlation is 
crucial to describe first-order nonmagnetic Mott transi- 
tions. The D-H attractive projector, if properly chosen, 
considerably reduces the variational energy for U W. 

(2) We clarified the optimized weight of the D-H bind- 
ing factor /aM, which rapidly decreases for r ^ 3, but 
becomes almost constant for r ^ 3 in the metallic regime. 
Thereby, it becomes clear why the simple conventional 
short-range D-H factors capture the essence of the Mott 
transition. In the insulating regime, the D-H binding fac- 
tor becomes very short-ranged; consequently, charge fluc- 
tuation, i.e. D-H pairs, is confined within r 2. 

(3) The D-D repulsive factors improve the energy only 
slightly, especially in the metallic regime, and do not in- 
duce a Mott transition by itself. However, the improve- 
ment in the insulating regime contributes to some down- 
ward shift of the Mott critical point U c /t. 

(4) Motivated by the Mott transition in the completely 
D-H bound state, we have renewed the picture of Mott 
transitions. Two characteristic length scales <dh and I dd 
are introduced; Idh broadly represents the size of a D-H 
pair, and ^dd the minimum distance between two dou- 
blons. The two lengths generally depend largely on U/t, 
and should be appropriately estimated. The Mott crit- 
ical point determined by the condition Idh = ^dd is 
consistent with the values estimated from other quanti- 
ties. This picture is applicable to a wide range of Mott 
transitions including the Bose Hubbard model. 39 - ) 

We leave some intriguing subjects for future studies: 
(i) Improvement of the critical point U c /t. (ii) Introduc- 
tion of explicit AF and superconducting correlation into 
the one-body part. 56 ) (hi) How the new picture works 
for doped cases, namely, doped Mott insulators like the 
cuprate superconductors. 
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